home *** CD-ROM | disk | FTP | other *** search
/ Turnbull China Bikeride / Turnbull China Bikeride - Disc 2.iso / STUTTGART / LANG / ADA / GNAT / !gcc / adainc / 4 / adb / s-arit64 < prev    next >
Text File  |  1996-02-12  |  18KB  |  678 lines
  1. ------------------------------------------------------------------------------
  2. --                                                                          --
  3. --                         GNAT COMPILER COMPONENTS                         --
  4. --                                                                          --
  5. --                      S Y S T E M . A R I T H _ 6 4                       --
  6. --                                                                          --
  7. --                                 B o d y                                  --
  8. --                                                                          --
  9. --                            $Revision: 1.5 $                              --
  10. --                                                                          --
  11. -- GNAT is free software;  you can  redistribute it  and/or modify it under --
  12. -- terms of the  GNU General Public License as published  by the Free Soft- --
  13. -- ware  Foundation;  either version 2,  or (at your option) any later ver- --
  14. -- sion.  GNAT is distributed in the hope that it will be useful, but WITH- --
  15. -- OUT ANY WARRANTY;  without even the  implied warranty of MERCHANTABILITY --
  16. -- or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License --
  17. -- for  more details.  You should have  received  a copy of the GNU General --
  18. -- Public License  distributed with GNAT;  see file COPYING.  If not, write --
  19. -- to  the Free Software Foundation,  59 Temple Place - Suite 330,  Boston, --
  20. -- MA 02111-1307, USA.                                                      --
  21. --                                                                          --
  22. -- As a special exception,  if other files  instantiate  generics from this --
  23. -- unit, or you link  this unit with other files  to produce an executable, --
  24. -- this  unit  does not  by itself cause  the resulting  executable  to  be --
  25. -- covered  by the  GNU  General  Public  License.  This exception does not --
  26. -- however invalidate  any other reasons why  the executable file  might be --
  27. -- covered by the  GNU Public License.                                      --
  28. --                                                                          --
  29. -- GNAT was originally developed  by the GNAT team at  New York University. --
  30. -- It is now maintained by Ada Core Technologies Inc (http://www.gnat.com). --
  31. --                                                                          --
  32. ------------------------------------------------------------------------------
  33.  
  34.  
  35.  
  36. with Interfaces; use Interfaces;
  37. with Unchecked_Conversion;
  38.  
  39. package body System.Arith_64 is
  40.  
  41.    pragma Suppress (Overflow_Check);
  42.    pragma Suppress (Range_Check);
  43.  
  44.    subtype Uns64 is Unsigned_64;
  45.    function To_Uns is new Unchecked_Conversion (Int64, Uns64);
  46.    function To_Int is new Unchecked_Conversion (Uns64, Int64);
  47.  
  48.    subtype Uns32 is Unsigned_32;
  49.  
  50.    -----------------------
  51.    -- Local Subprograms --
  52.    -----------------------
  53.  
  54.    function "+" (A, B : Uns32) return Uns64;
  55.    function "+" (A : Uns64; B : Uns32) return Uns64;
  56.    pragma Inline ("+");
  57.    --  Length doubling additions
  58.  
  59.    function "-" (A : Uns64; B : Uns32) return Uns64;
  60.    pragma Inline ("-");
  61.    --  Length doubling subtraction
  62.  
  63.    function "*" (A, B : Uns32) return Uns64;
  64.    function "*" (A : Uns64; B : Uns32) return Uns64;
  65.    pragma Inline ("*");
  66.    --  Length doubling multiplications
  67.  
  68.    function "/" (A : Uns64; B : Uns32) return Uns64;
  69.    pragma Inline ("/");
  70.    --  Length doubling division
  71.  
  72.    function "rem" (A : Uns64; B : Uns32) return Uns64;
  73.    pragma Inline ("rem");
  74.    --  Length doubling remainder
  75.  
  76.    function "&" (Hi, Lo : Uns32) return Uns64;
  77.    pragma Inline ("&");
  78.    --  Concatenate hi, lo values to form 64-bit result
  79.  
  80.    function Lo (A : Uns64) return Uns32;
  81.    pragma Inline (Lo);
  82.    --  Low order half of 64-bit value
  83.  
  84.    function Hi (A : Uns64) return Uns32;
  85.    pragma Inline (Hi);
  86.    --  High order half of 64 bit value
  87.  
  88.    ---------
  89.    -- "+" --
  90.    ---------
  91.  
  92.    function "+" (A, B : Uns32) return Uns64 is
  93.    begin
  94.       return Uns64 (A) + Uns64 (B);
  95.    end "+";
  96.  
  97.    function "+" (A : Uns64; B : Uns32) return Uns64 is
  98.    begin
  99.       return A + Uns64 (B);
  100.    end "+";
  101.  
  102.    ---------
  103.    -- "-" --
  104.    ---------
  105.  
  106.    function "-" (A : Uns64; B : Uns32) return Uns64 is
  107.    begin
  108.       return A - Uns64 (B);
  109.    end "-";
  110.  
  111.    ---------
  112.    -- "*" --
  113.    ---------
  114.  
  115.    function "*" (A, B : Uns32) return Uns64 is
  116.    begin
  117.       return Uns64 (A) * Uns64 (B);
  118.    end "*";
  119.  
  120.    function "*" (A : Uns64; B : Uns32) return Uns64 is
  121.    begin
  122.       return A * Uns64 (B);
  123.    end "*";
  124.  
  125.    ---------
  126.    -- "/" --
  127.    ---------
  128.  
  129.    function "/" (A : Uns64; B : Uns32) return Uns64 is
  130.    begin
  131.       return A / Uns64 (B);
  132.    end "/";
  133.  
  134.    -----------
  135.    -- "rem" --
  136.    -----------
  137.  
  138.    function "rem" (A : Uns64; B : Uns32) return Uns64 is
  139.    begin
  140.       return A rem Uns64 (B);
  141.    end "rem";
  142.  
  143.    ---------
  144.    -- "&" --
  145.    ---------
  146.  
  147.    function "&" (Hi, Lo : Uns32) return Uns64 is
  148.    begin
  149.       return Shift_Left (Uns64 (Hi), 32) or Uns64 (Lo);
  150.    end "&";
  151.  
  152.    --------
  153.    -- Hi --
  154.    --------
  155.  
  156.    function Hi (A : Uns64) return Uns32 is
  157.    begin
  158.       return Uns32 (Shift_Right (A, 32));
  159.    end Hi;
  160.  
  161.    --------
  162.    -- Lo --
  163.    --------
  164.  
  165.    function Lo (A : Uns64) return Uns32 is
  166.    begin
  167.       return Uns32 (A and 16#FFFF_FFFF#);
  168.    end Lo;
  169.  
  170.    --------------------------
  171.    -- Add_With_Ovflo_Check --
  172.    --------------------------
  173.  
  174.    function Add_With_Ovflo_Check (X, Y : Int64) return Int64 is
  175.       R : constant Int64 := To_Int (To_Uns (X) + To_Uns (Y));
  176.  
  177.    begin
  178.       if X >= 0 then
  179.          if Y < 0 or else R >= 0 then
  180.             return R;
  181.          end if;
  182.  
  183.       else -- X < 0
  184.          if Y > 0 or else R < 0 then
  185.             return R;
  186.          end if;
  187.       end if;
  188.  
  189.       raise Constraint_Error;
  190.    end Add_With_Ovflo_Check;
  191.  
  192.    -----------------------------
  193.    -- Divide_With_Ovflo_Check --
  194.    -----------------------------
  195.  
  196.    function Divide_With_Ovflo_Check (X, Y : Int64) return Int64 is
  197.    begin
  198.       if Y = 0
  199.         or else (Y = (-1) and then X = (-(2 ** 63)))
  200.       then
  201.          raise Constraint_Error;
  202.       else
  203.          return X / Y;
  204.       end if;
  205.    end Divide_With_Ovflo_Check;
  206.  
  207.    -------------------
  208.    -- Double_Divide --
  209.    -------------------
  210.  
  211.    procedure Double_Divide
  212.      (X, Y, Z : Int64;
  213.       Q, R    : out Int64;
  214.       Round   : Boolean)
  215.    is
  216.       Xu  : constant Uns64 := To_Uns (abs X);
  217.       Xhi : constant Uns32 := Hi (Xu);
  218.       Xlo : constant Uns32 := Lo (Xu);
  219.  
  220.       Yu  : constant Uns64 := To_Uns (abs Y);
  221.       Yhi : constant Uns32 := Hi (Yu);
  222.       Ylo : constant Uns32 := Lo (Yu);
  223.  
  224.       Zu  : constant Uns64 := To_Uns (abs Z);
  225.       Zhi : constant Uns32 := Hi (Zu);
  226.       Zlo : constant Uns32 := Lo (Zu);
  227.  
  228.       T1, T2     : Uns64;
  229.       Du, Qu, Ru : Uns64;
  230.       Den_Pos    : Boolean;
  231.  
  232.    begin
  233.       if Yu = 0 or else Zu = 0 then
  234.          raise Constraint_Error;
  235.       end if;
  236.  
  237.       --  Compute Y * Z. Note that if the result overflows 64 bits unsigned,
  238.       --  then the rounded result is clearly zero (since the dividend is at
  239.       --  most 2**63 - 1, the extra bit of precision is nice here!)
  240.  
  241.       if Yhi /= 0 then
  242.          if Zhi /= 0 then
  243.             Q := 0;
  244.             R := X;
  245.             return;
  246.          else
  247.             T2 := Yhi * Zlo;
  248.          end if;
  249.  
  250.       else
  251.          if Zhi /= 0 then
  252.             T2 := Ylo * Zhi;
  253.          else
  254.             T2 := 0;
  255.          end if;
  256.       end if;
  257.  
  258.       T1 := Ylo * Zlo;
  259.       T2 := T2 + Hi (T1);
  260.  
  261.       if Hi (T2) /= 0 then
  262.          Q := 0;
  263.          R := X;
  264.          return;
  265.       end if;
  266.  
  267.       Du := Lo (T2) & Lo (T1);
  268.       Qu := Xu / Du;
  269.       Ru := Xu rem Du;
  270.  
  271.       --  Deal with rounding case
  272.  
  273.       if Round and then Ru > (Du - 1) / 2 then
  274.          Qu := Qu + 1;
  275.       end if;
  276.  
  277.       --  Set final signs (RM 4.5.5(27-30))
  278.  
  279.       Den_Pos := (Y < 0) = (Z < 0);
  280.  
  281.       --  Case of dividend (X) sign positive
  282.  
  283.       if X >= 0 then
  284.          R := To_Int (Ru);
  285.  
  286.          if Den_Pos then
  287.             Q := To_Int (Qu);
  288.          else
  289.             Q := -To_Int (Qu);
  290.          end if;
  291.  
  292.       --  Case of dividend (X) sign negative
  293.  
  294.       else
  295.          R := -To_Int (Ru);
  296.  
  297.          if Den_Pos then
  298.             Q := -To_Int (Qu);
  299.          else
  300.             Q := To_Int (Qu);
  301.          end if;
  302.       end if;
  303.    end Double_Divide;
  304.  
  305.    -------------------------------
  306.    -- Multiply_With_Ovflo_Check --
  307.    -------------------------------
  308.  
  309.    function Multiply_With_Ovflo_Check (X, Y : Int64) return Int64 is
  310.       Xu  : constant Uns64 := To_Uns (abs X);
  311.       Xhi : constant Uns32 := Hi (Xu);
  312.       Xlo : constant Uns32 := Lo (Xu);
  313.  
  314.       Yu  : constant Uns64 := To_Uns (abs Y);
  315.       Yhi : constant Uns32 := Hi (Yu);
  316.       Ylo : constant Uns32 := Lo (Yu);
  317.  
  318.       T1, T2 : Uns64;
  319.  
  320.    begin
  321.       if Xhi /= 0 then
  322.          if Yhi /= 0 then
  323.             raise Constraint_Error;
  324.          else
  325.             T2 := Xhi * Ylo;
  326.          end if;
  327.  
  328.       else
  329.          if Yhi /= 0 then
  330.             T2 := Xlo * Yhi;
  331.          else
  332.             return X * Y;
  333.          end if;
  334.       end if;
  335.  
  336.       T1 := Xlo * Ylo;
  337.       T2 := T2 + Hi (T1);
  338.  
  339.       if Hi (T2) /= 0 then
  340.          raise Constraint_Error;
  341.       end if;
  342.  
  343.       T2 := Lo (T2) & Lo (T1);
  344.  
  345.       if X >= 0 then
  346.          if Y >= 0 then
  347.             return To_Int (T2);
  348.          else
  349.             return -To_Int (T2);
  350.          end if;
  351.       else -- X < 0
  352.          if Y < 0 then
  353.             return To_Int (T2);
  354.          else
  355.             return -To_Int (T2);
  356.          end if;
  357.       end if;
  358.  
  359.    end Multiply_With_Ovflo_Check;
  360.  
  361.    -------------------
  362.    -- Scaled_Divide --
  363.    -------------------
  364.  
  365.    procedure Scaled_Divide
  366.      (X, Y, Z : Int64;
  367.       Q, R    : out Int64;
  368.       Round   : Boolean)
  369.    is
  370.       Xu  : constant Uns64 := To_Uns (abs X);
  371.       Xhi : constant Uns32 := Hi (Xu);
  372.       Xlo : constant Uns32 := Lo (Xu);
  373.  
  374.       Yu  : constant Uns64 := To_Uns (abs Y);
  375.       Yhi : constant Uns32 := Hi (Yu);
  376.       Ylo : constant Uns32 := Lo (Yu);
  377.  
  378.       Zu  : Uns64 := To_Uns (abs Z);
  379.       Zhi : Uns32 := Hi (Zu);
  380.       Zlo : Uns32 := Lo (Zu);
  381.  
  382.       D1, D2, D3, D4 : Uns32;
  383.       --  The dividend, four digits (D1 is high order)
  384.  
  385.       Q1, Q2 : Uns32;
  386.       --  The quotient, two digits (Q1 is high order)
  387.  
  388.       S1, S2, S3 : Uns32;
  389.       --  Value to subtract, three digits (S1 is high order)
  390.  
  391.       Qu : Uns64;
  392.       Ru : Uns64;
  393.       --  Unsigned quotient and remainder
  394.  
  395.       Scale : Natural;
  396.       --  Scaling factor used for multiple-precision divide. Dividend and
  397.       --  Divisor are multiplied by 2 ** Scale, and the final remainder
  398.       --  is divided by the scaling factor. The reason for this scaling
  399.       --  is to allow more accurate estimation of quotient digits.
  400.  
  401.       T1, T2, T3 : Uns64;
  402.       --  Temporary values
  403.  
  404.    begin
  405.       --  First do the multiplication, giving the four digit dividend
  406.  
  407.       T1 := Xlo * Ylo;
  408.       D4 := Lo (T1);
  409.       D3 := Hi (T1);
  410.  
  411.       if Yhi /= 0 then
  412.          T1 := Xlo * Yhi;
  413.          T2 := D3 + Lo (T1);
  414.          D3 := Lo (T2);
  415.          D2 := Hi (T1) + Hi (T2);
  416.  
  417.          if Xhi /= 0 then
  418.             T1 := Xhi * Ylo;
  419.             T2 := D3 + Lo (T1);
  420.             D3 := Lo (T2);
  421.             T3 := D2 + Hi (T1);
  422.             T3 := T3 + Hi (T2);
  423.             D2 := Lo (T3);
  424.             D1 := Hi (T3);
  425.  
  426.             T1 := (D1 & D2) + (Xhi * Yhi);
  427.             D1 := Hi (T1);
  428.             D2 := Lo (T1);
  429.  
  430.          else
  431.             D1 := 0;
  432.          end if;
  433.  
  434.       else
  435.          if Xhi /= 0 then
  436.             T1 := Xhi * Ylo;
  437.             T2 := D3 + Lo (T1);
  438.             D3 := Lo (T2);
  439.             D2 := Hi (T1) + Hi (T2);
  440.  
  441.          else
  442.             D2 := 0;
  443.          end if;
  444.  
  445.          D1 := 0;
  446.       end if;
  447.  
  448.       --  Now it is time for the dreaded multiple precision division. First
  449.       --  an easy case, check for the simple case of a one digit divisor.
  450.  
  451.       if Zhi = 0 then
  452.          if D1 /= 0 or else D2 >= Zlo then
  453.             raise Constraint_Error;
  454.  
  455.          --  Here we are dividing at most three digits by one digit
  456.  
  457.          else
  458.             T1 := D2 & D3;
  459.             T2 := Lo (T1 rem Zlo) & D4;
  460.  
  461.             Qu := Lo (T1 / Zlo) & Lo (T2 / Zlo);
  462.             Ru := T2 rem Zlo;
  463.          end if;
  464.  
  465.       --  If divisor is double digit and too large, raise error
  466.  
  467.       elsif (D1 & D2) >= Zu then
  468.          raise Constraint_Error;
  469.  
  470.       --  This is the complex case where we definitely have a double digit
  471.       --  divisor and a dividend of at least three digits. We use the classical
  472.       --  multiple division algorithm (see  section (4.3.1) of Knuth's "The Art
  473.       --  of Computer Programming", Vol. 2 for a description (algorithm D).
  474.  
  475.       else
  476.          --  First normalize the divisor so that it has the leading bit on.
  477.          --  We do this by finding the appropriate left shift amount.
  478.  
  479.          Scale := 0;
  480.  
  481.          if (Zhi and 16#FFFF0000#) = 0 then
  482.             Scale := 16;
  483.             Zu := Shift_Left (Zu, 16);
  484.          end if;
  485.  
  486.          if (Hi (Zu) and 16#FF00_0000#) = 0 then
  487.             Scale := Scale + 8;
  488.             Zu := Shift_Left (Zu, 8);
  489.          end if;
  490.  
  491.          if (Hi (Zu) and 16#F000_0000#) = 0 then
  492.             Scale := Scale + 4;
  493.             Zu := Shift_Left (Zu, 4);
  494.          end if;
  495.  
  496.          if (Hi (Zu) and 16#C000_0000#) = 0 then
  497.             Scale := Scale + 2;
  498.             Zu := Shift_Left (Zu, 2);
  499.          end if;
  500.  
  501.          if (Hi (Zu) and 16#8000_0000#) = 0 then
  502.             Scale := Scale + 1;
  503.             Zu := Shift_Left (Zu, 1);
  504.          end if;
  505.  
  506.          Zhi := Hi (Zu);
  507.          Zlo := Lo (Zu);
  508.  
  509.          --  Note that when we scale up the dividend, it still fits in four
  510.          --  digits, since we already tested for overflow, and scaling does
  511.          --  not change the invariant that (D1 & D2) >= Zu.
  512.  
  513.          T1 := Shift_Left (D1 & D2, Scale);
  514.          D1 := Hi (T1);
  515.          T2 := Shift_Left (0 & D3, Scale);
  516.          D2 := Lo (T1) or Hi (T2);
  517.          T3 := Shift_Left (0 & D4, Scale);
  518.          D3 := Lo (T2) or Hi (T3);
  519.          D4 := Lo (T3);
  520.  
  521.          --  Compute first quotient digit. We have to divide three digits by
  522.          --  two digits, and we estimate the quotient by dividing the leading
  523.          --  two digits by the leading digit. Given the scaling we did above
  524.          --  which ensured the first bit of the divisor is set, this gives an
  525.          --  estimate of the quotient that is at most two too high.
  526.  
  527.          if D1 = Zhi then
  528.             Q1 := 2 ** 32 - 1;
  529.          else
  530.             Q1 := Lo ((D1 & D2) / Zhi);
  531.          end if;
  532.  
  533.          --  Compute amount to subtract
  534.  
  535.          T1 := Q1 * Zlo;
  536.          T2 := Q1 * Zhi;
  537.          S3 := Lo (T1);
  538.          T1 := Hi (T1) + Lo (T2);
  539.          S2 := Lo (T1);
  540.          S1 := Hi (T1) + Hi (T2);
  541.  
  542.          --  Adjust quotient digit if it was too high
  543.  
  544.          loop
  545.             exit when S1 < D1;
  546.  
  547.             if S1 = D1 then
  548.                exit when S2 < D2;
  549.  
  550.                if S2 = D2 then
  551.                   exit when S3 <= D3;
  552.                end if;
  553.             end if;
  554.  
  555.             Q1 := Q1 - 1;
  556.  
  557.             T1 := (S2 & S3) - Zlo;
  558.             S3 := Lo (T1);
  559.             T1 := (S1 & S2) - Zhi;
  560.             S2 := Lo (T1);
  561.             S1 := Hi (T1);
  562.          end loop;
  563.  
  564.          --  Subtract from dividend (note: do not bother to set D1 to
  565.          --  zero, since it is no longer needed in the calculation).
  566.  
  567.          T1 := (D2 & D3) - S3;
  568.          D3 := Lo (T1);
  569.          T1 := (D1 & Hi (T1)) - S2;
  570.          D2 := Lo (T1);
  571.  
  572.          --  Compute second quotient digit in same manner
  573.  
  574.          if D2 = Zhi then
  575.             Q2 := 2 ** 32 - 1;
  576.          else
  577.             Q2 := Lo ((D2 & D3) / Zhi);
  578.          end if;
  579.  
  580.          T1 := Q2 * Zlo;
  581.          T2 := Q2 * Zhi;
  582.          S3 := Lo (T1);
  583.          T1 := Hi (T1) + Lo (T2);
  584.          S2 := Lo (T1);
  585.          S1 := Hi (T1) + Hi (T2);
  586.  
  587.          loop
  588.             exit when S1 < D2;
  589.  
  590.             if S1 = D2 then
  591.                exit when S2 < D3;
  592.  
  593.                if S2 = D3 then
  594.                   exit when S3 <= D4;
  595.                end if;
  596.             end if;
  597.  
  598.             Q2 := Q2 - 1;
  599.  
  600.             T1 := (S2 & S3) - Zlo;
  601.             S3 := Lo (T1);
  602.             T1 := (S1 & S2) - Zhi;
  603.             S2 := Lo (T1);
  604.             S1 := Hi (T1);
  605.          end loop;
  606.  
  607.          T1 := (D3 & D4) - S3;
  608.          D4 := Lo (T1);
  609.          T1 := (D2 & Hi (T1)) - S2;
  610.          D3 := Lo (T1);
  611.  
  612.          --  The two quotient digits are now set, and the remainder of the
  613.          --  scaled division is in (D3 & D4). To get the remainder for the
  614.          --  original unscaled division, we rescale this dividend.
  615.  
  616.          Qu := Q1 & Q2;
  617.          Ru := Shift_Right (D3 & D4, Scale);
  618.       end if;
  619.  
  620.       --  Deal with rounding case
  621.  
  622.       if Round and then Ru > (Zu - 1) / 2 then
  623.          Qu := Qu + 1;
  624.       end if;
  625.  
  626.       --  Set final signs (RM 4.5.5(27-30))
  627.  
  628.       --  Case of dividend (X * Y) sign positive
  629.  
  630.       if (X >= 0 and then Y >= 0)
  631.         or else (X < 0 and then Y < 0)
  632.       then
  633.          R := To_Int (Ru);
  634.  
  635.          if Z > 0 then
  636.             Q := To_Int (Qu);
  637.          else
  638.             Q := -To_Int (Qu);
  639.          end if;
  640.  
  641.       --  Case of dividend (X * Y) sign negative
  642.  
  643.       else
  644.          R := -To_Int (Ru);
  645.  
  646.          if Z > 0 then
  647.             Q := -To_Int (Qu);
  648.          else
  649.             Q := To_Int (Qu);
  650.          end if;
  651.       end if;
  652.  
  653.    end Scaled_Divide;
  654.  
  655.    -------------------------------
  656.    -- Subtract_With_Ovflo_Check --
  657.    -------------------------------
  658.  
  659.    function Subtract_With_Ovflo_Check (X, Y : Int64) return Int64 is
  660.       R : constant Int64 := To_Int (To_Uns (X) - To_Uns (Y));
  661.  
  662.    begin
  663.       if X >= 0 then
  664.          if Y > 0 or else R >= 0 then
  665.             return R;
  666.          end if;
  667.  
  668.       else -- X < 0
  669.          if Y <= 0 or else R < 0 then
  670.             return R;
  671.          end if;
  672.       end if;
  673.  
  674.       raise Constraint_Error;
  675.    end Subtract_With_Ovflo_Check;
  676.  
  677. end System.Arith_64;
  678.